In this notebook we will build maps for the species of Anchylorhynchus, using Natural Earth as base map and combining species with compatible distributions in a single map.
rm(list=ls())
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.1 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.3 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(rnaturalearth)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
library(RStoolbox)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.0.1
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
Let’s read specimen data, keep only georreferrenced data and remove duplicates:
sp_data = readxl::read_excel('especimes estudados_2020_rev2.xlsx')
sp_data = sp_data %>%
filter(str_detect(updated_locality,'mistake',negate=T) | is.na(updated_locality)) %>%
dplyr::select(accepted_name,latitude, longitude) %>%
filter(!is.na(latitude)&!is.na(longitude)) %>%
distinct %>%
arrange(accepted_name, latitude, longitude) %>%
rename(species = accepted_name)
sp_data
Now, let’s download the base map (commented because only run once).
#basemap = ne_download(scale = 'large', type = 'GRAY_HR_SR_OB_DR', category = 'raster', destdir = './map')
After map is downloaded, let’s crop it to the region of interest:
basemap = raster::stack('map/GRAY_HR_SR_OB_DR/GRAY_HR_SR_OB_DR.tif')
bbox = raster::extent(c(xmin=-95,xmax=-25,ymin=-45,ymax=19))
cropped = raster::crop(basemap,bbox)
map_layer = ggRGB(cropped, ggLayer = TRUE,r = 1,g=1,b=1)
#the following is just to visualize the layer
print(ggRGB(cropped, r = 1,g=1,b=1))
countries = ne_countries(scale=50, returnclass = 'sf') %>%
st_transform(crs = 4326) %>%
st_crop(c(xmin = bbox@xmin,
xmax = bbox@xmax,
ymin = bbox@ymin,
ymax = bbox@ymax))
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
Now, let’s load a table indicating which species should be plotted together.
sp_map = readr::read_csv('map_combination.csv')
## Parsed with column specification:
## cols(
## species = col_character(),
## map = col_double()
## )
sp_map
Now, let’s create a function that plots each map:
plot_map = function(spp){
this_occ = sp_data %>%
filter(species %in% spp$species) %>%
mutate(species = str_c('A. ',species)) %>%
st_as_sf(coords=c('longitude', 'latitude'), crs=4326)
ggplot(this_occ) +
map_layer +
geom_sf(color='black',fill=NA, size=0.2, data=countries) +
geom_sf(aes(shape = species, fill = species), size = 1.5, show.legend = "point") +
coord_sf(expand = F, )+
scale_fill_brewer(type='qual') +
scale_shape_manual(values = 21:24) +
#guides(colour = guide_legend(override.aes = list(size=0.25))) +
theme(
axis.text = element_text(size=7),
axis.ticks = element_line(),
#legend.text = element_text(face='italic',size=6,margin=margin(t=.1,r=0,b=.1,l=.1,unit='pt')),
legend.text = element_text(face='italic',size=6),
legend.margin = margin(t=0,r=2,b=2,l=2,unit='pt'),
legend.justification = c(1,1),
legend.position = c(1,1),
legend.title= element_blank(),
legend.box.margin = margin(t=0,r=0,b=0,l=0,unit='pt'),
legend.background = element_rect(fill = 'white',colour = NA),
legend.box.background = element_rect(fill = 'white',colour = NA),
legend.key = element_rect(fill = 'white',colour = NA),
legend.key.height = unit(8,'pt'),
legend.key.width = unit(5,'pt')
) +
xlim(bbox@xmin,bbox@xmax) +
ylim(bbox@ymin,bbox@ymax)
}
Now let’s plot all maps:
plots = split(sp_map,sp_map$map) %>%
purrr::map(plot_map)
plots = 1:length(plots) %>%
purrr::map(~plots[[.x]] + labs(tag=LETTERS[.x]))
plots
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
And now arrange them in a page:
p = grid.arrange(arrangeGrob(grobs=plots,
nrow = 3))
ggsave('maps_2.pdf', plot=p, device='pdf',width=168,height=200,units='mm',useDingbats=F)
ggsave('maps_2.eps', plot=p, device='eps',width=168,height=200,units='mm')
plot(p)